# function that sets up the model and runs the counterfactual
function set_up_model(
    ramp,
    amenities,
    produc,
    realloc_workers,
    realloc_pollution,
    realloc_production,
    abatement,
    congestion,
    opt,
    vsl_elastic,
    info,
    heterog,
    tol,
    θ = 4.,
    α = .274, 
    γ = .481,
    ι = 1.,
    ξ_mult = 1.
    )
    """
    wrapper function that
    sets up and runs 
    the entire model
    """

    ########################
    ### set up parameters
    ########################

    # read in pollutant and productivity response
    pollutant_response = -1 .* Matrix(CSV.read("data/counterfactual/inputs/emissions_reductions.csv", DataFrame))[1:5]
    # NOx/PM25/SO2/VOC elasticities are the manufacturing industry median from Shapiro Walker 2018 AER
    # NH3 is not in their paper so we assign it to the median of PM2.5 which is a catch-all for small particulates
    ξ = (nh3 = .0023 .* ξ_mult, nox = .0038 .* ξ_mult, pm25 = .0023 .* ξ_mult, so2 = .0028 .* ξ_mult, voc = .0068 .* ξ_mult)

    if congestion
        ζ_c = 0.3
        ζ_a = 0.2
    else
        ζ_c = 0.0
        ζ_a = 0.0
    end

    if vsl_elastic
        ϵ_vsl = 0.4
    else
        ϵ_vsl = 0.0
    end

    p = (
        cf_year = 1990,    # counterfactual year
        year = 1997,       # year
        γ = γ,         # aggregate labor share: https://www.bls.gov/opub/mlr/2017/article/estimating-the-us-labor-share.htm
        θ = θ,            # trade elasticity (Simonovska and Waugh 2014)
        α = α,          # US manufacturing consumption share (Rudik et al 2021)
        ι = ι,            # migration elasticity (Jaworski et al)
        ζ_c = ζ_c,             # congestion elasticity (Allen and Arkolakis 2013)
        ζ_a = ζ_a,             # agglomeration elasticity (Allen and Arkolakis 2013)
        ξ = [ξ.nh3, ξ.nox, ξ.pm25, ξ.so2, ξ.voc],  # pollution elasticity 
        ϵ_vsl = ϵ_vsl,          # income elasticity of VSL
        vsl_inc = 67.521,        # base income for VSL calculation in $10,000
        damp = ι == 1 ? .1 : 0.25,          # damping
        tol = tol,        # convergence tolerance
        β_T = produc, # effect on productivity
        β_E = pollutant_response,       # change in emission price
        ramp = ramp, # ramp up nonattainment
        amenities = amenities,  # amenities impacts (Bhat update)
        realloc_workers = realloc_workers,     # mobility adjustments (μ update)
        realloc_pollution = realloc_pollution, # pollution transport (ap3 model)
        realloc_production = realloc_production, # changing market potential (ma update)
        abatement = abatement, # emission price effects (η change in counterfactual)
        congestion = congestion, # congestion
        opt = opt, # optimal run
        info = info, # households have information about amenities,
        heterog = heterog, # heterogenous treatment effect of nonattainment (deprecated)
        ξ_mult = ξ_mult, # multiplier on pollution elasticities
        run_base = true # run base model
        ) 

    # load data
    data = load_data(p)

    # solve model
    @time solve_counterfactual(p, data)

end